clear all;
clc;

% ----------------------- Geometry Parameters -----------------------------

% import data
patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

patnum = 5;

epidur  = 1; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = 0; % 0, -10, -20

% location of electrode

eletag_epi = 'epi_';
fdist_elec_epi = cordgeom(patnum,8);

eletag_sub = 'sub_';
if(subloc==1)
    fdist_elec_sub = cordgeom(patnum,9);
elseif(subloc==2)
    fdist_elec_sub = 1-cordgeom(patnum,9);
else
    return;
end

ptag=['p',num2str(patnum),'_'];
loctag_epi = ['d',strrep(num2str(fdist_elec_epi),'.','p')];
loctag_sub = ['d',strrep(num2str(fdist_elec_sub),'.','p')];

% ----------------------- Spine Geometry ----------------------------------

% cord geometry
delr_dura = 0.2;

betageom = patientgeom(patnum,:);
a_cord = betageom(1)/2;
b_cord = betageom(2)/2;
a_csf  = betageom(3)/2-delr_dura;
b_csf  = betageom(4)/2-delr_dura;
a_dura = a_csf + delr_dura;
b_dura = b_csf + delr_dura;
ycenter = 13.61;
csf_center = [0,ycenter];
ytopepi = 25;

y_cord = (ycenter-b_csf) + betageom(5)*2*b_csf;
cord_center = [0,y_cord];

% ------------------------ Electrode Location -----------------------------

ytmp = ycenter+b_csf+delr_dura;
yelec_epi = ytmp+(ytopepi-ytmp)*fdist_elec_epi;

yelec_sub = (y_cord+b_cord)+(ycenter+b_csf-(y_cord+b_cord))*fdist_elec_sub;

elec_center_epi = [0,yelec_epi];
elec_center_sub = [0,yelec_sub];

th_rad_0 = 0*(pi/180);
th_rad_n10 = -10*(pi/180);
th_rad_n20 = -20*(pi/180);
Arot_0   = [cos(th_rad_0),-sin(th_rad_0);sin(th_rad_0),cos(th_rad_0)];
Arot_n10 = [cos(th_rad_n10),-sin(th_rad_n10);sin(th_rad_n10),cos(th_rad_n10)];
Arot_n20 = [cos(th_rad_n20),-sin(th_rad_n20);sin(th_rad_n20),cos(th_rad_n20)];

% epidural
rp_epi_0 = Arot_0*[0;yelec_epi-y_cord];
xe_epi_0 = rp_epi_0(1); ye_epi_0 = rp_epi_0(2) + y_cord;
relec_epi_0 = [xe_epi_0,ye_epi_0];

rp_epi_n10 = Arot_n10*[0;yelec_epi-y_cord];
xe_epi_n10 = rp_epi_n10(1); ye_epi_n10 = rp_epi_n10(2) + y_cord;
relec_epi_n10 = [xe_epi_n10,ye_epi_n10];

rp_epi_n20 = Arot_n20*[0;yelec_epi-y_cord];
xe_epi_n20 = rp_epi_n20(1); ye_epi_n20 = rp_epi_n20(2) + y_cord;
relec_epi_n20 = [xe_epi_n20,ye_epi_n20];

% subdural
rp_sub_0 = Arot_0*[0;yelec_sub-y_cord];
xe_sub_0 = rp_sub_0(1); ye_sub_0 = rp_sub_0(2) + y_cord;
relec_sub_0 = [xe_sub_0,ye_sub_0];

rp_sub_n10 = Arot_n10*[0;yelec_sub-y_cord];
xe_sub_n10 = rp_sub_n10(1); ye_sub_n10 = rp_sub_n10(2) + y_cord;
relec_sub_n10 = [xe_sub_n10,ye_sub_n10];

rp_sub_n20 = Arot_n20*[0;yelec_sub-y_cord];
xe_sub_n20 = rp_sub_n20(1); ye_sub_n20 = rp_sub_n20(2) + y_cord;
relec_sub_n20 = [xe_sub_n20,ye_sub_n20];


% white matter
th_unit = 0:pi/20:2*pi;
rcord   = a_cord*b_cord./sqrt((b_cord*cos(th_unit)).^2+(a_cord*sin(th_unit)).^2); 
xcord   = cord_center(1) + rcord.*cos(th_unit);
ycord   = cord_center(2) + rcord.*sin(th_unit);
% csf space
rcsf   = a_csf*b_csf./sqrt((b_csf*cos(th_unit)).^2+(a_csf*sin(th_unit)).^2); 
xcsf   = csf_center(1) + rcsf.*cos(th_unit);
ycsf   = csf_center(2) + rcsf.*sin(th_unit);
% dura
rdura   = a_dura*b_dura./sqrt((b_dura*cos(th_unit)).^2+(a_dura*sin(th_unit)).^2); 
xdura   = csf_center(1) + rdura.*cos(th_unit);
ydura   = csf_center(2) + rdura.*sin(th_unit);

% grey matter
x1 = [0.00 ,0.62 ,1.24 ,1.87 ,2.09 ,2.24 ,1.45 ,2.03 ,2.22 ,2.03 ,1.21 ,0.62 ,0.00];
y1 = [-0.44,-0.12,0.96,2.04,2.06,1.88,0.43,-1.12,-1.62,-2.12,-2.18,-1.61,-0.88];
y1 = y1+cord_center(2);
xgrey = [x1,-fliplr(x1)];
ygrey = [y1,fliplr(y1)];
         
% ----------------------- Importing Data ----------------------------------

pol = -1;
noaxons = 200;

prefix_epi = [ptag,'neg_',eletag_epi,'bp_'];
prefix_sub = [ptag,'neg_',eletag_sub,'bp_'];

cd(['patient',num2str(patnum),'_thresholds']);

% DC fibers
dcthresh_epi_0 = load([prefix_epi,'dc_act_',loctag_epi,'_t0','.txt']);
[dcthresh_epi_0(:,1),dc_si_epi_0] = sort( dcthresh_epi_0(:,1)*pol );
dcthresh_sub_0 = load([prefix_sub,'dc_act_',loctag_sub,'_t0','.txt']);
[dcthresh_sub_0(:,1),dc_si_sub_0] = sort( dcthresh_sub_0(:,1)*pol );

dcthresh_epi_n10 = load([prefix_epi,'dc_act_',loctag_epi,'_tn10','.txt']);
[dcthresh_epi_n10(:,1),dc_si_epi_n10] = sort( dcthresh_epi_n10(:,1)*pol );
dcthresh_sub_n10 = load([prefix_sub,'dc_act_',loctag_sub,'_tn10','.txt']);
[dcthresh_sub_n10(:,1),dc_si_sub_n10] = sort( dcthresh_sub_n10(:,1)*pol );

dcthresh_epi_n20 = load([prefix_epi,'dc_act_',loctag_epi,'_tn20','.txt']);
[dcthresh_epi_n20(:,1),dc_si_epi_n20] = sort( dcthresh_epi_n20(:,1)*pol );
dcthresh_sub_n20 = load([prefix_sub,'dc_act_',loctag_sub,'_tn20','.txt']);
[dcthresh_sub_n20(:,1),dc_si_sub_n20] = sort( dcthresh_sub_n20(:,1)*pol );

dc_xy = load([ptag,'dc_xyloc.txt']);
dc_xy = dc_xy + ( cord_center'*ones(1,noaxons) )';

dc_xy_epi_0 = dc_xy(dc_si_epi_0,:);
dc_xy_epi_n10 = dc_xy(dc_si_epi_n10,:);
dc_xy_epi_n20 = dc_xy(dc_si_epi_n20,:);
dc_xy_sub_0 = dc_xy(dc_si_sub_0,:);
dc_xy_sub_n10 = dc_xy(dc_si_sub_n10,:);
dc_xy_sub_n20 = dc_xy(dc_si_sub_n20,:);

% DR fibers
drthresh_epi_0 = load([prefix_epi,'dr_act_',loctag_epi,'_t0','.txt']);
[drthresh_epi_0(:,1),dr_si_epi_0] = sort( drthresh_epi_0(:,1)*pol );
drthresh_sub_0 = load([prefix_sub,'dr_act_',loctag_sub,'_t0','.txt']);
[drthresh_sub_0(:,1),dr_si_sub_0] = sort( drthresh_sub_0(:,1)*pol );

drthresh_epi_n10 = load([prefix_epi,'dr_act_',loctag_epi,'_tn10','.txt']);
[drthresh_epi_n10(:,1),dr_si_epi_n10] = sort( drthresh_epi_n10(:,1)*pol );
drthresh_sub_n10 = load([prefix_sub,'dr_act_',loctag_sub,'_tn10','.txt']);
[drthresh_sub_n10(:,1),dr_si_sub_n10] = sort( drthresh_sub_n10(:,1)*pol );

drthresh_epi_n20 = load([prefix_epi,'dr_act_',loctag_epi,'_tn20','.txt']);
[drthresh_epi_n20(:,1),dr_si_epi_n20] = sort( drthresh_epi_n20(:,1)*pol );
drthresh_sub_n20 = load([prefix_sub,'dr_act_',loctag_sub,'_tn20','.txt']);
[drthresh_sub_n20(:,1),dr_si_sub_n20] = sort( drthresh_sub_n20(:,1)*pol );

cd ..;

% ----------------------- DC Activated before DR --------------------------

nocase = 6;
xy_selact = cell(1,nocase);
relec = zeros(nocase,2);

nodc_epi_0   = sum( dcthresh_epi_0(:,1) < drthresh_epi_0(1,1) );
nodc_epi_n10 = sum( dcthresh_epi_n10(:,1) < drthresh_epi_n10(1,1) );
nodc_epi_n20 = sum( dcthresh_epi_n20(:,1) < drthresh_epi_n20(1,1) );
nodc_sub_0   = sum( dcthresh_sub_0(:,1) < drthresh_sub_0(1,1) );
nodc_sub_n10 = sum( dcthresh_sub_n10(:,1) < drthresh_sub_n10(1,1) );
nodc_sub_n20 = sum( dcthresh_sub_n20(:,1) < drthresh_sub_n20(1,1) );

xy_selact{1} = dc_xy_epi_0(1:nodc_epi_0,:);
xy_selact{2} = dc_xy_epi_n10(1:nodc_epi_n10,:);
xy_selact{3} = dc_xy_epi_n20(1:nodc_epi_n20,:);
xy_selact{4} = dc_xy_sub_0(1:nodc_sub_0,:);
xy_selact{5} = dc_xy_sub_n10(1:nodc_sub_n10,:);
xy_selact{6} = dc_xy_sub_n20(1:nodc_sub_n20,:);

relec(1,:) = relec_epi_0;
relec(2,:) = relec_epi_n10;
relec(3,:) = relec_epi_n20;
relec(4,:) = relec_sub_0;
relec(5,:) = relec_sub_n10;
relec(6,:) = relec_sub_n20;

% ----------------------- Spatial Activation ------------------------------

figure;
for k = 1:nocase
    
    perdc = 100*length(xy_selact{k})/noaxons;

    subplot(2,3,k);
    hold on;
    fill(xcord,ycord,[1,1,1],'LineWidth',3);
    plot((xcsf+xdura)/2,(ycsf+ydura)/2,'k--','LineWidth',4);
    fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
    plot(relec(k,1),relec(k,2),'ko','MarkerSize',10,...
         'MarkerFaceColor','w','LineWidth',4);
    plot(xy_selact{k}(:,1),xy_selact{k}(:,2),'k.','MarkerSize',25);
    set(gca,'XTick',[],'YTick',[]);
    

    ttl = [num2str(perdc,2),' %'];
    title(ttl,'FontSize',36);
    axis tight;
    axis square;
    hold off;
    
end

figure;
subplot(2,3,k);
hold on;
fill(xcord,ycord,[1,1,1],'LineWidth',3);
plot((xcsf+xdura)/2,(ycsf+ydura)/2,'k--','LineWidth',4);
fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
plot(relec(k,1),relec(k,2),'ko','MarkerSize',12,...
     'MarkerFaceColor','w','LineWidth',4);
plot(xy_selact{k}(:,1),xy_selact{k}(:,2),'k.','MarkerSize',25);
set(gca,'XTick',[],'YTick',[]);
axis off;

L1 = 'white matter'; 
L2 = 'dura mater'; 
L3 = 'grey matter'; 
L4 = 'electrode location';
L5 = 'DC fibers activated';
legend(L1,L2,L3,L4,L5,'Location','BO');
set(gca,'FontSize',36);